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ABSTRACT 

It  is  shown  that  by  using  the  simplest  construction  of  discrete  dipoles, 

the  operation  count  for  solving  the  Dirichlet  problem  of  Poisson's  equation 

2 

by  the  capacitance  matrix  method  does  not  exceed  constant  times  n log  n, 
n = 1/h  for  certain  first  and  second  order  schemes  of  interpolating  boundary 
conditions. 


AMS(MOS)  Subject  Classification:  65N20 

Key  Words:  Capacitance  matrix,  Green’s  function, 

Poisson  equation,  conjugate  gradient  method, 
collectively  conpact  operator. 
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SIGNIFICANCE  AND  EXPLANATION 

The  Dirichlet  problem  for  the  Poisson  equation  is  the  following:  Given  a 

function  f and  a function  g , find  a fiinction  u such  that 

u + u = f on  fi  , 

XX  yy 

u = g on  3 n. 

Here  is  a simply  connected  domain  with  boundary  3 R . 

The  problem  has  wide  applications  in  electrostatics,  elasticity,  teit5)er- 
ature  distributions  and  plasma  physics.  Its  solution  by  finite  difference  or 
finite  elements  methods  have  received  considerable  attention.  It  is 
known  that  if  R is  a rectangle,  then  fast  Fourier  transform  methods  are  very 
efficient  in  solving  the  linear  system  of  equations  arising  from  finite  differ- 
ence or  finite  element  discretizations. 

There  seems  to  be  no  such  short  cut  to  the  solution  of  these  equations 

when  R is  a general  region.  In  many  conventional  methods,  the  operation 

3/2 

count  is  usually  proportional  to  N (N  is  the  number  of  mesh  points  in 

R ) while  at  least  N computer  storage  is  required.  These  methods  are 
therefore  undesirable  when  N is  very  large.  In  this  paper  we  describe  an 
algorithm  and  prove  mathematically  that  the  operation  count  of  this  algorithm 
can  be  proportional  to  N log  N.  While  some  versions  of  our  algorithm  also 
require  at  least  N computer  storage , there  is  one  version  that  requires  less 
them  N/3  computer  storage. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summeury  lies  with  MRC,  and  not  with  the  author  of  this  report. 


FAST  POISSON  SOLVERS  ON  GENERAL  TWO  DIMENSIONAL  REGIONS 
FOR  THE  DIRICHLET  PROBLEM 
A.  S.  L.  Shieh 

§1 . Introduction 

Over  the  past  ten  years,  very  fast  nun>erical  methods  have  been  developed  to  solve 
Poisson's  or  Helmholtz's  equation  on  certain  simple  regions  with  Dirichlet,  Neumann  or 
periodic  boundary  conditions.  See  e.g.  12],  [3],  [8],  [9],  [12],  [19]  and  [21].  These 
methods  can  only  be  used  for  regions  and  boundary  conditions  that  allow  for  separation  of 
the  variables.  Typical  examples  are  Poisson's  or  Helmholtz's  equations  in  Cartesian  co- 
ordinates on  rectangular  regions  with  boundary  conditions  that  do  not  change  type  along  any 
of  the  sides  of  the  rectangle.  In  these  special  cases,  the  operation  count  for  solving  the 
discrete  problem  is  almost  proportional  to  the  number  of  mesh  points. 

The  purpose  of  this  paper  is  to  establish  similar  results  for  the  Poisson  equation 
on  general  regions.  In  this  work  we  are  only  concerned  with  finite  difference  schemes 
of  first  and  second  order  accuracy  for  the  Dirichlet  problem  on  simply  connected  bounded 
domains  with  smooth  boundaries.  A formal  discrete  potential  theory  motivated  by  the  clas- 
sical potential  theory  is  incorporated  into  the  so-called  capacitance  matrix  method.  It  is 

shown  that  by  using  the  simplest  construction  of  discrete  dipoles  in  our  Ansatz,  it  is 

2 

possible  to  have  an  algorithm  the  operation  count  of  which  is  proportional  to  N log  N , 
where  h = 1/n  is  the  mesh  size.  Some  numerical  results  are  given  in  section  9 and  a 
brief  survey  of  past  work  in  this  direction  is  given  in  section  8. 


Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-75-C-0024  and  the  Energy 
Research  eind  Development  Administration. 


§2 . Certain  results  from  classical  potential  theory. 

We  give  only  a very  brief  review  of  a few  results  of  classical  potential  theory.  For  a 

detailed  exposition  see  e.g.  (10,,  [14J  , [25]  and  [28].  We  define  the  potential  / resulting 

from  a charge  distribution  r on  a smooth  boundary  curve  3(2  by 

(x)  = (lA)/  p(C)  logiids(C)  . 

2 2 2 

Here  x = (x^^.x^),  C=  ” “ ^''l”^l'  * Green's  function 

(1/2T')  log  T which  we  shall  denote  by  G*  satisfies 
A(1/2ti)  log  r = 6{x)  , 

where  6(x)  is  the  delta  function.  Similarly  the  potential  W of  a dipole  density  h on 
3(2  defined  by 

(2.1)  '•(x)  = (l/Ti)/  U(C)  3g*/3v_  ds(C)  . 

3(2  ^ 

We  adopt  here  the  convention  that  the  normal  direction  of  3(2  is  towards  the  exterior  of 

the  region  (2  in  which  we  wemt  to  solve  our  problem. 

The  interior  Dirichlet  problem  can  be  reduced  to  a Fredholm  integral  equation  of  the 

second  kind  if  we  make  the  double  layer  Ansatz  as  follows.  Let 

u(x)  = -(1/211)  //  f(5)  logr  dC  + d/n)  / h (C)  3GV3Vr  ds(£)  = u (x)  +lv'(x)  , 

(2  3(2  ® 

for  the  solution  of 

-Au  = f,  X e (2 

(2.2) 

u = g,  X e 3(2 

It  can  be  shown  that  the  dipole  density  U satisfy  the  following  integral  equation 

(2.3)  p + (1/11)  / pOG/3(2j.)ds  = g - u I = g . 

3(2  ® 3(2 

This  is  a well  posed  problem  of  the  form 

(2.4)  (I  + K)p  « g , 

where  K is  a compact  operator  defined  by  the  integral  above. 

If  we  instead  atten^it  to  use  a single  layer  Ansatz  for  the  Dirichlet  problem  we  obtain  a 
Fredholm  integral  equation  of  the  first  kind.  It  has  the  form 

/(x)  ” g ■ “g  I an  ' ' 

which  is  an  ill  posed  problem. 

To  illustrate  the  distribution  of  the  eigenvalues  of  the  conjiact  operator  K 
in  equation  (2.4),  we  study  the  case  when  fi  is  an  ellipse  with 
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Y = (a-b)/(a+b)  where  a and  b are  the  half  axes  of  the  ellipse.  It  is  known  (see  e.g. 
p.  135  of  [25])  that  K = and 

(2.5)  X^(K+k'^)  = , i - 1,2 

On  the  other  hand,  both  the  interior  and  exterior  Neumann  problems  can  also  be  reduced  to 
Fredholm  integral  equations  of  the  second  kind  if  we  make  the  single  layer  Ansatz.  The 

charge  density  p for  the  exterior  Neumann  problem  satisfies 

T * 

(I  + K )p  = g , 

* 

for  some  suitably  chosen  function  g defined  on  3f).  The  existence  and  uniqueness  problems 
for  the  solution  of  equation  (2.3)  can  therefore  be  determined  from  that  of  equation  (2.4) 

h 

and  vice  versa.  Finally,  we  remark  that  the  G in  equations  (2.1)  and  (2.3)  can  be  re- 
placed by  the  Green's  function  on  a sufficiently  large  rectangle  with  zero  Dirichlet  bound- 
ary conditions. 
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§3.  The  capacitance  matrix  method 


In  this  section  we  develop  a similar,  formal  potential  theory  for  the  discrete 
problems  arising  from  the  original  Dirichlet  problem  (2.2).  See  also  Sections  3 and  4 of  [29  1 
for  a similar  discussion.  We  shall  assume  that  uniform  mesh  sizes  in  both  coordinate  directions 
are  used. 

We  replace  the  Laplace  operator  by  the  five-point  formula.  The  fundamental  solu- 
tion (1/211)  log(l/r),  used  in  Section  2,  will  be  replaced  by  its  discrete  analogue,  the  dis- 
crete Green's  function  on  the  entire  plane,  which  we  shall  denote  by  G.  Properties  and 
efficient  methods  of  generating  G and  its  undivided  differences  will  be  studied  in  Section  4. 
An  efficient  method  of  computing  Gv  for  arbitrary  N N vectors  v is  also  given 

in  Section  4.  We  will  denote  by  B the  matrix  representing  the  five-point  discrete  Laplacian 
2 

h using  undivided  differences,  on  the  entire  plane.  We  then  divide  the  set  of  mesh  points 

into  three  disjoint  sets  12,  ,3fi.  and  (.CU)^.  The  set  9f2,  contains  all  the  irregular  mesh 

n h h h 

points  in  S2,  i.e.  mesh  points  that  do  not  have  all  four  neighlaours  within  the  open  set  Q . 

Q.  is  the  set  of  regular  mesh  points  inside  U and  (Cfl),  contains  the  remaining,  the  exte- 
h h 

rior  mesh  points. 

We  then  set  up  the  matrix  equation 
(3.1)  Au  = V 

that  we  are  solving  as  follows.  We  use  the  same  discretization  formula  for  Ijoth  A and  B 
on  U (cn)j^.  For  points  in  a linear  combination  of  the  discrete  Laplacian  and  inter- 

polation formulas  of  first  or  second  order  accuracy  for  the  boundary  conditions  are  used.  The 
values  of  the  solution  at  the  exterior  mesh  points  are  always  eliminated  from  the  discrete 
Laplacicin,  centered  at  an  irregular  mesh  point.  This  guarantees  that  A is  a reducible  matrix 
with  no  couplings  to  the  exterior  mesh  points  from  the  irregular  mesh  points.  If  P is  a 
suitably  chosen  permutation  matrix,  then 

T 

PAP 

where  A is  the  coefficient  matrix  for  our  discrete  problem  on  O.  u 
11  h 

seen  that  the  solution  on  12.  U 312.  will  not  be  influenced  by  either  the 

h h 

data  on  (Cl2)  . 

h 
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3(2jj.  It  is  easily 
solution  or  the 

I 


I 


The  matrix  A differs  from  B by  only  m rows  where  m is  the  cardinal  number  of 

. We  can  therefore  write 
h 

T 

A = B + U W , 

where  the  matrices  U and  W have  m columns.  The  matrix  U represents  an  extension 

operator.  It  maps  any  mesh  function  defined  only  on  into  a function  on  all  mesh  points. 

T 

Its  transpose,  U , is  a trace  operator  mapping  any  mesh  function  defined  for  all  mesh  points 
into  its  restriction  to  We  easily  verify  that 


W 


T 


= u'^(A-B)  . 


We  now  describe  our  method  for  solving  the  discrete  problem  (3.1).  Guided  by  the  contin- 
uous analog  we  ma)ce  the  Ansatz, 

(3.2)  u = Gv  + GV  P . 

The  vector  Gv  satisfies  BGv  = v.  The  m-vector  p is  determined  by  solving  a system  of 
mxm  linear  equations  derived  below.  The  mesh  function  VP  should  vanish  on  Each 

column  of  the  matrix  V represents  a discrete  dipole  of  unit  strength.  Let  such  a column 

T 

corresponding  to  P e 3(1^^  be  regarded  as  a mesh  function,  denoted  by  V (P)  . We  require  that 
T 

(l/h^)V  (p)  u = [3u/3v] (P)  + 0(h).  Here  h^  = h/cosu  where  a is  the  angle  between  the  normal 
through  the  irregular  mesh  point  P and  the  closest  coordinate  eixis.  In  particular,  if 
that  the  western  and  northwestern  neighbours  of  P in  (CS))^^,  then 


(3.3)  [V^u] (P)  = u(P)  - (1  - tan  a)  u (w)  - (tan  a)u(NW). 

We  now  use  our  Ansatz  and  compute  the  residual  vector, 

(3.4)  Au-v  = (B  + uw'^)  (GF  + GV)J)  - Fv 

= (V  + uw'^GVjy  + Uw'^GFv. 

From  the  properties  of  U and  V,  it  follows  that  the  residuals  are  zero  for  all  x € 

T 

To  derive  a linear  system  of  equations  for  the  vector  P we  multiply  equation  (3.4)  by  U . 

T T 

It  is  easy  to  verify  that  U U = I and  U v = I . Here  I is  the  raxm  identity  matrix. 

m mm 

We  thus  obtain 

(3.5)  (I  + W^GV)p  = -w'^GF  . 

m 

This  choice  of  p ma)<e  the  residuals  zero  for  all  x e Hence  substitution  of  V*  in 

equation  (3.2)  will  provide  us  a solution  on  Li  351  if  equation  (3.5)  is  solvable.  Note 

h h 

that  the  residuals  will  in  general  not  be  equal  to  zero  for  all  x f matrix  on 
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the  lett-hand  side  of  Equation  (3.5)  is  the  capacitance  matrix  C . We  shall  refer  to 

Equation  (3.5)  as  the  capacitance  ’^•atrix  equation. 

T -1 

In  the  special  case  when  v = u U v,  we  can  simply  make  the  Ansatz  u = B vy.  It  is 
easily  seen  that  the  residual  Au-v  will  again  be  zero  at  x e The  capacitance  matrix 

equation  now  becomes 

(3.6)  Cy  = u’’v  . 

If  Equation  (3.6)  is  solvable,  then  Au  = v will  also  be  zero  on  3 il  . The  solvability  of 

h 

Equations  (3.5)  and  (3.6)  will  be  discussed  in  Section  6. 

We  now  describe  our  choices  of  difference  equations  at  the  irregular  mesh  points.  We 

approximate  the  boundary  conditions  by  interpolation  schemes  of  first  or  second  order  accuracy, 

which  we  shall  refer  to  as  schemes  I?,  Ib  and  II  respectively. 

We  start  with  Scheme  II.  Let  P e and  P*  be  its  closest  point  on  Bfl.  Let 

W,E,N  and  S l>e  the  western,  eastern,  northern  and  southern  neighbours  of  P on  the  mesh 

respectively.  We  assume  that  the  local  orientation  of  the  boundary  is  such  that  either  both 

W and  N are  in  (CQ),  or  only  W is  in  (Cf2).  . Assume  that  both  W and  N are  in  (Cfl),  . 

h h h 

Let  dj^  denote  h^^/h  where  hj^  is  the  distance,  along  a mesh  line  parallel  with  the  x^-axis, 
between  the  mesh  point  P and  the  boundary  312.  Hence  dj^  e (0,1).  The  Dirichlet  data  at 
this  point  on  3fl  is  denoted  by  u . The  values  of  d_  and  u are  similarly  defined.  We 

then  approximate  u^^  and  u^^  by  (1/2)  [ (1+d^) u (W)  + (l-dj^)u(E)]  and  (1/2)  (l+d2)u(N)  + 
(l-d2)u(S)l  respectively.  By  combining  the  above  with  the  five-point  formula  for  the  Laplacian 
and  eliminating  u(W)  and  u(N)  between  them,  we  obtain 

(3.6)  4u(P)  - :2dj/(l+d^)]u(E)  - (2  d^/d+dj" 

= h^f(P)  + [2/(l+d,)]u,  + [2/(l+d  )]u  . 

1 W z N 

If  only  W is  in  then  we  obtain 

(3.7)  4u(P)  - [2  dj^/(l+dj^)]u(E)  - u(N)  - u(S) 

- h^f  (P)  + [2/(l+d^)]u„  . 

We  now  describe  the  two  variants  of  Scheme  I,  namely  Scheme  la  and  Ib.  In  Scheme  Ib, 

if  both  W emd  N are  in  (Cfi),  , we  obtain 

h 

' 4u(P)  - u(S)  - u(E)  = h^f(P)  + u,  + u„  . 

W N 
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If  only  W is  in  (C;. ),  , we  obtain 
h 

(3.9)  2[4u(P)  - u(S)  - 11(E)  - u(N)l  = 2 (h^f(P)  + . 

The  scalinq  factor  2 is  largely  artifical  and  is  put  in  only  for  the  convenience  of 

theoretical  estimates  in  Sections  5 and  6.  In  Scheme  la,  we  seek  to  eliminate  this  scaling 

factor  while  retaining  the  theoretical  convenience.  We  require  that  Equation  (3.8)  should  be 

used  regardless  of  whether  both  W and  N are  in  (CQ)^  or  only  W is  in  (CQ),  . The 

h h 

matrices  for  all  the  above  three  schemes  are  of  positive  type.  Hence,  the  results  in 

{4]  or  [13]  apply  and  all  these  schemes  are  convergent. 

There  is  an  important  alternative  to  the  above  approach.  Instead  of  the  discrete  Green's 
function  of  the  entire  plane,  we  may  use  the  discrete  analog  of  the  Green's  function  on  a 

sufficiently  large  square  S with  zero  boundary  conditions  as  our  G in  equations  (3.2)  and 

-1  T 

(3.5).  In  this  case  G = ; A = + UW  . Here  denotes  the  matrix  representing  the 

2 

discrete  Laplacian  h A.  on  S and  zero  boundary  values  on  the  grid  points  of  3S.  The 
h 

T -1  . 

residual  Au-v  will  again  be  zero  on  UudQ  if  C = U A B V is  nonsinqular, 

h h u 

Finally  we  come  to  the  central  question  as  to  whether  the  capacitance  matrix  equation 

(3.5)  is  closely  related  to  the  Fredholm  integral  equation  (2.2)?  It  is  Itnown  (see  e.g.  [16]) 

that  the  conjugate  gradient  method  converges  superl inearly  for  Fredholm  integral  equations  of 

the  second  )cind.  In  our  experiments  we  normally  fail  to  observe  superlinear  convergence.  To 

understand  this  fully,  we  split  up  the  matrices  C into  two  parts  as  follows. 

C = B + K,  . 
h h 


The  matrices  are  defined  by 


(3.10) 


Bj^(P,Q) 


C(P,Q),  if  d(P,Q)  < /h 


= 0 , otherwise  . 

They  are  therefore  the  near  diagonal  parts  of  C ; and  the  matrices  1^  are  the  remaining 
parts,  the  off  diagonal  parts  of  C . 

It  will  be  shown  in  Section  6 that  for  Schemes  La  and  Ib,  and,  after  a suitable  scaling, 

for  schemes  II,  the  matrices  K,  are  closely  related  to  the  compact  integral  operator  K in 

h 

Equations  (2.2)  or  (2.3).  The  matrices  Bj^,  however,  will  not  in  general  be  formal  approxi- 
mations to  the  identy  operator.  In  fact,  the  algebraic  row  sums  of  B^^  need  not  always  )?e 
equal  to  one. 
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It  is»  however,  shown  in  (16)  that  for  operator  equations  with  symmetric  positive  defi- 


nite operators  of  the  form  B + K with  B symmetric  positive  definite  and  K symmetric 
compact,  the  conjugate  gradient  method  will  converge  linearly  with  asymptotic  rate  of  con- 
vergence governed  only  by  the  spectral  condx  .ion  number  of  B . We  therefore  proceed  to 

study  the  special  condition  numbers  of  B.  in  Section  5 and  the  distribution  of  singular 

h 

values  of  K,  and  C in  Section  6.  We  shall  show  in  Section  7 that  the  asymptotic  con- 
h 

vergence  of  the  conjugate  gradient  method  for  solving  the  capacitance  matrix  equations  will 

depend  essentially  on  the  spectral  condition  number  of  B,  . 

h 

We  now  discuss  briefly  two  different  methods  of  irrplementing  our  algorithm  and  the 

operation  count  involved.  We  use  the  conjugate  gradient  method  to  solve 

T T 

C Cu  = C b, 

where  b denotes  the  right  hand  side  of  the  capacitance  matrix  equation.  The  solution  u 
is  then  confuted  from  (3.2).  In  the  first  method,  we  generate  the  capacitance  matrix 
explicitly.  Assume  that  the  G in  equations  (3.2)  and  (3.5)  is  the  discrete  Green's  func- 
tion on  the  entire  plane.  Because  of  translational  invariance  it  is  only  necessary  to  com- 
pute G with  the  second  parameter  fixed  at  the  origin.  It  is  shown  in  section  4 that  only 
one  call  of  fast  Poisson  solved  on  a sufficiently  large  rectangle  is  needed  to  generate  G 

and  only  two  calls  of  a similar  solver  is  needed  to  compute  the  final  solution  and  the  right 

2 

hand  side  b . The  operation  count  of  the  algorithm  is  therefore  constant  N log  N 
2 

+ 2 C^m  + 0(m),  where  is  the  number  of  iterations  needed  to  achieve  a certain  accuracy. 

If  the  G in  equation  (3.5)  is  B^^  , it  is  desirable  to  use  the  second  method  where 
the  solution  \i  is  computed  by  an  iterative  implicit  method  first  appearing  in  (151.  The 
operation  count  for  confuting  y is  proportional  to  CQ(iTH'm^)N  , where  m^  is  the  number  of 
nonzero  entries  in  the  matrix  V provided  that  a special  fast  solver  is  used  in  the  process. 
See  section  4 of  [311  for  details.  It  will  be  shown  in  section  7 that  cannot  exceed 

constant  log  m if  G = B^^^  is  used  in  (3.5)  for  all  domains  with  sufficiently  smooth  bound- 
aries and  that  is  uniformly  bounded  in  some  special  cases  if  the  discrete  Green's  func- 

tion on  the  entire  plane  is  used  in  (3.5).  The  total  operation  count  of  our  algorithm  there- 

2 

fore  does  not  exceed  constant  N log  N . 
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§4.  Properties  and  fast  generation  of  G , the  discrete  Green's  function  on  the  entire  plane , 

T 

and  the  efficient  computations  of  Gv,  GVn  and  W Gv. 

A discrete  fundamental  solution  of  the  five-point  Laplacian  wiuh  respect  to  the 

origin  is  a mesh  function  y that  satisfies 

r if  U = 0 , 

(4.1)  \Y(yh)  = J 

1.  0 if  p 0 , 

where  U has  integer  components  and  p^  . 

Clearly  Y is  unique  up  to  an  arbitrary  linear  function.  The  constants  involved  will  be 
chosen  so  that  we  have  a proper  discrete  analog  of  the  logarithmic  potential.  The  resulting 
discrete  fundamental  solution  will  then  be  our  discrete  Green's  function  G. 

It  is  established  in  [27]  that  if  g(r,s)  denotes 

(4.2)  (2/11)  /^[l-cos(sA)  exp(- [r]  p)  ]/sinhp  dl 

0 

where 

(4.3)  cos  A + cosh  p = 2 

(4.4)  p/A  >-1  as  p 0, 
then  the  function  G defined  by 

(4.5)  G(r)c,sh)  = (1/4)  g(r,s)  + (1/271)109  h - (1/471)  (log  8 + 2y^) 

is  the  desired  Green's  function  of  the  entire  plane.  Here  Yj^  is  the  Euler's  constant;  r 
and  s integers  and  h is  the  mesh  size. 

It  is  shown  in  [27]  that 

(4.6)  g(0,0)  = 0,  g(0,l)  = 1 

(4.7)  g(r,s)  = g(s,r)  = g(-s,r)  = g(s,-r) , 

(4.8)  g(r,r)  = (4/71)  [1  + 1/3  + ...  + l/(2r-l)l  , 

(4.9)  g(r,s)  = (I/ti) log (s^+r^)  + (l/7i)  (log  8 + 2Yj)  + o(l/r), 

(4.10)  g(r,s)  - g(r,t)  = (l/ii)  log  ( (r^+s^) /(t^+r^)  1 + od/r). 

We  have  found  it  necessary  to  obtain  sharper  estimates  for  the  remainder  term  in  (4.10) 
when  t = s+1  and  a similar  estimate  for  g(r+l,s)  - g(r,s) . 

Theorem  4.1.  Let  r,s  and  t be  nonnegative  integers  with  r 8;  s » t-1.  Then 
(4.11a)  G(rh,sh)  - G(rh,th)  = (1/477) log [ (s^+r^)/ ( t^+r^)  ) + Rj^(r,s)  + R2(r,s)  + Rj(r,s), 

where 
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(4.11b)  R,(r,s)  = (1/2411)  / (X  -rX'^)sin  [ (s  + 1/2)  ' 1 

^ 0 


(4.11c)  R (r,s)  = -(1/211)  J {[(l/30)>i'*+(7/96)r'^+(l/288)r^X®]  . e'^^^sin  ( (2s+l ) V2]  , 

0 2/3 

(4. lid)  |R2(r,s)|  < (0.'>7)r‘^e”^'^^'^  + (167)r'^  . 

Proof.  We  hiave 

(ii,'2)  ig (r ,s) -g (r , t)  ] = / lcos(tX)  - cos(sX)]e  '^'^[sinhy]  ^ dX  , 

0 

and  it  is  Icnovm  that 

/“’[cos(tX)  - cos(sX)]e”’^^X"^dX  = (1/2)  log(  (s^+r^)/ (t^+r^)  ] . 


Since 


(4.12)  / |cos(tX)  - cos(sX)|  e '^^X  ^dX  < (2/it)r 

IT 

it  suffices  to  estimate 

/^[cos(sX)  - cos(tA))[e  "^^X  ^ - e '^'^(sinhy)  ^]dX. 

0 

The  integrand  in  the  above  expression  will  be  denoted  by  J.  Let  e = (1.5)r  , r > 8.  We 

have 

(4.13)  f^J  d X = tcos(sX)-cos(tX)  ]e"’^^[2  sin  (X/2)  ] "^F  (X)  dX  , 

0 0 

where 

(4.14)  F(X)  S 2[sin(X/2)]X'^-2  sin(X/2)  [sinhpl’^e'^"^'’^  ■ 

By  (4.3)  and  (4.4),  we  have 

(4.15)  sinh  y = 2Isin(X/2)]  (l+sin^(X/2)]^'^^  , 

(4.16)  e"^  = 2 - cosX  -2[sin(X/2)]  (1  + sin^  (X/2)  ] . 

It  is  easily  verified  that  for  0 £ X <^  £, 

(4.17)  = 1 + (1/12)  (X^+X^)  + (1/288)X®  + C^(X)X’, 

with  |Cj^(X)|  £ 0.035,  0 £,X^1.  Hence, 

e(X-y)r  ^ (1/12)  (rX^  + rX^)  + (l/288)r^X^  + C2(r,X)  . 

Here  |C2(r,X)|  £ (0.04)rX^  + (O.Ol)r^X^  + (O.OOl)r^X^  . Therefore, 

F(X)  = (1/12)  [X^-rX^  - (7/8)rX®]  - (1/30)X‘*  - (l/288)r^X^  + Cj(r,X). 

Here  lCj(r,X)|  £ (O.Ol)X®  + (0.051)rX^  + (0.0104)r^X®  + (O.OOl)r^X^  . 

Hence, 

(4.18)  /^'JdX  = 2ii[Rj^(r,s)  + R2(r,s))  + Zj^(r,s)+  Z2(r,s), 

0 

where  (4.19) 

(4.19)  |z  (r,s)|  < r e"’^^C  (r,X)dX 

0 

< 1047  r-’^ 
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(4.20) 


lz^(r,s)  1 < (l/6)r'^ 


2/? 


It  remains 

to  estimate  J^j  dA  . 

Clearly , 

t 

(4.21) 

1 

1 J JdA 

! ^ / “tA  -rii . 

1 ^ J (e  + e )d>-  , 

It  can  be 

c 

shown  that 

f 

“P 

e us  a decreasing  function  of  A for 

!33l.  It 

is  also  shown  there  that 

(4.22) 

e-^  < 

-0.76 

e , X = 

0.8 

< 

-0.91 

e , X = 

1 

< 

-1.07  , 

e , X = 

1.2 

< 

-1.31  , 

e , X = 

11/2. 

By  (4.17), 

e-^  = e 

r^(l  + C*(X)X)  , ]c* 

(X) 1 <0.1,  0 < X < 0 

so  that 

e'""  1 

e-°-"'"  , 0 < X < 

0.8. 

Hence, 

,11  -rp 

J e 

dX  < (0.32)e”°'^^’^ 

0.8 

and 

(4.23) 

-rA 

+ e JdA  < (2.15)r 

-1  -1.35r^'^^  -1  -1 

e -re 

By  combining  (4.12),  (4.19),  (4.20)  and  (4.22),  we  see  that 

, ,«:  2/3 

(4.24)  211  iRjCr.s)!  j<  2 . 3 r e + 1047  r 

The  theorem  then  follows  from  (4.5),  (4.13),  (4.18)- (4.21)  and  (4 . 23- (4 . 24) . 

Theorem  4.2.  Let  r,s  and  t be  nonnegative  integers  with  t ^ 8;  t = r+1 . Then 
(4.25a) 


G(rh,sh)  - G(th,sh)  = (1/4TI)  log(  (r^+s^)/ (t^+s^)  1 + S^(r,s)  + S2(r,s)  + S^(r,s), 


where 

(4.25b) 

(4.25c) 


(4.25d) 


S (r,s)  = /“‘(1/24TT)  [l^-rX^+X^e'^^^]e"‘*"*^^^^'^cos  sX  (sX) 

0 

S2(r,s)  = - /"  (1/211)  ( (l/30)X^+(7/96)rX^  + (l/288)r^X^  + (1/12)  X'^e'^'^^l  • 

° -(r+l/2)X 

• e cos  sA(sA) 


|Sj(r,s)|  < (1.9)r"®  + (206)r'’  + (0.5)r'^ 


2/3 
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Proof . 


We  have 


(7^/2)  (g(r  ,s)-g(t,s)  I = / [e  ^^-e  ^coE(sA)dA  + J j dA , 


_ , wt  *tM,  . , »-l  -tA.-l  -rU , . . -1  . -rA^-I, 

J - cos(sA)[e  (sinhp)  - e A - e (sinhvi)  + e A ) 


It  is  known  that 


(4.26)  P cos(sA)dA  = (1/2)  log  I (r^+s'^)  / (t'^+s'^)  ] . 

0 

Clearly, 

(4.27)  /“  < (l/n)  . 

-7T  * 

It  therefore  suffices  to  estimate  J J dX.  Let 
. * . 0 

(4.28)  J =^1+^2' 


(4.29)  = cos(sX)[e  ^-e  ’^'^(sinhy)  ^)  (1-e  ^)  , 

(4.30)  S cos(sX)[e  *^-e  ^]e  "^^X  ^ . 

Let  e=(1.5)r"^^^  . By  (4.16), 

(l-e‘^)  [2  sin(X/2]"^  = -sin(X/2)  + [1  + sin^ (X/2) ) < 1. 


-sin(X/2)  + (1  + sin^(X/2)  1^/^  = e~  + b^(X)X^  , 


|bj^(X)  1 < 1/12,  0 < X < e . 


Hence,  by  (4.14) 


/q  = /jj  cos(sX)e  ^^■'■^/2)X  ^ Ej^(r,s)  , 


E^(r,s)|  < /”  (1/12)X^  F(X)e"'^  dX  < (5.84)r"^  + 364r‘“  + lO^r 


8 . ..5  -10 


By  (4.17), 


(4.33)  / J - / ((1/12) (X+X)  + (1/288)X’  + C (X) X®] cos  (sX) e‘  dX, 

0^0  ^ 

with  |Cj^(X)|  ^ 0.035.  It  is  easily  seen  (see  also  p.  38  of  (301)  that 


I ,*  jil  ..  -1  -(1.5)r  2/3  ^ ,,,  -1  -1. 

IJ  J dX|  £ r e + (2.15)r  e 


-1.35r  2/3  -1  -rn 

- e . 


The  theorem  easily  follows  from  (4.19),  (4.20)  and  (4. 26- (4. 34) . 
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Theorems  4.1  and  4.2  provide  accurate  estimates  for  distant  values  of  the  undivided 
differences  of  G . One  can  use  those  estimates  as  Dirichlet  conditions  for  a fast  Poisson 
solver  to  generate  all  the  values  that  are  needed  to  set  up  the  capacitance  matrix  C . One 
can  also  use  integer  arithmetic  as  in  [27]  to  construct  a table  of  values  of  G for  r,s  < 7. 
See  e,g.  Table  II  on  p.  292  of  (27J  or  Table  I on  p.  41  of  [30] . By  using  the  above  tables 

and  Theorems  4.1  and  4.2,  we  obtain  the  following. 

Theorem  4.3.  Let  r and  t = s+1  be  positive  integers.  Then 

(4.35)  G(sh,rh)  - G(th,sh)  = (1/4TI) log [ (s^+r^) / (t^+r^)  ] + R(a,r) 


where 

(4.36)  (R(s,r)|  < (0.34)  min  {s  ^,r 


Moreover, 


7 7 7 

(4.37)  inax{  'i  max  |R(s,r)l,  \ max  |R(s,r)|,  1 max  [ R(s , r ) | } j<  0. 01 . 

r=2  s£r  s=2  r£s  s=r  r<s 

We  next  investigate  the  monotone  behavior  of  the  undivided  differences  of  G in  certain 
directions.  Let 


G^(i,j)  E G((i+l)h,jh)  - G(ih,jh) 

Gy(i,j)  E G(ih,(j+l)h)  -G(ih,jh)  . 

By  the  five-point  formula  and  symmetry, 

(4.38)  G (i-l,j)  = G (i,j-l)  i*0  or  j*0; 

XX  yy 

(4.39)  2G^(0,j)  = -G^(0,j-1)  jxO  . 

Theorem  4.4  Let  r and  s be  nonnegative  integers.  Then  G^(s,r),  G^is,r),  -G^^(s,r) 

and  G (l,r)  - G (0,r)  are  always  positive:  and  G (s-l,r)  is  always  nonnegative  for 

yy  yy  xx 

r ^ s, 

g£gsi-  Except  for  the  result  on  G^^(s-l,r),  the  proof  for  all  the  other  results  are  similar.  We 
first  estimate  the  values  of  the  expression  for  s = 0,  r > 8 using  Theorem  4.1  or  4.2.  By 
symmetry,  the  results  hold  for  s ^ 8,  r = 0 as  well.  We  then  verify  with  the  aid  of  Table  II 
on  p.  292  of  [271  that  the  same  results  hold  for  s = 0,  r 3 and  3^8,  r = 0.  Since  the 
five  point  formula  is  satisfied  at  all  points  rh  and  sh  with  r>n  and  s > 0,  an  application 
of  discrete  maximum  principle  immediately  yields  the  desired  result. 

The  proof  for  G^^(s-l,r)  is  as  follows.  By  symmetry  and  (4.38),  we  note  that 
G (s,r)  5 0 for  s » r . 

XX 
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By  symmetry. 


G (s-l,r)  « 2G  (0,r)  for  s ^ 0 . 

XX  X 

An  application  of  the  discrete  maximum  principle  therefore  completes  the  proof  of  the  theorem. 

This  concludes  our  discvssion  for  undivided  differences  of  G . We  now  proceed  to 
obtain  better  estimates  for  the  remainder  terms  in  (4.8)  and  (4.9). 

The  following  theorem  is  an  immediate  consequence  of  (4.8)  and  some  well  known  result  of 


asymptotic  series.  See  also  p.  325  of  (5  ]. 

Theorem  4.5  Let  r be  any  positive  integer.  Then 
(4.40)  G(rh,rh)  = {1/4tt) log (2r^h^)  + (l/48Ti)r'^  + r’  , 

where 


(4.41) 


< {7/1920ir)r 


-4 


Theorem  4.6  Let  r > s;  r and  s are  nonnegative  integers.  Then 

(4.42)  G(rh,sh)  = (l/4Tr)  log  f (s^+r^)  h^J  - (1/2411)  (r^+s^)  + (1/311)  r^s^  (r^+s^)  + L(r,s 

^ere 

2/3 

(4.43)  |L(r,s)|  < (2.51/ir)r"^  + (154/ll)r”^+  (l/il)log  r 

Proof.  As  in  the  proof  of  Theorem  4.1,  we  have 

/ J dX  = /^[cOs(sX)  - cos(tX))e  *^^(2  sin(X/2)l  ^F(X)dX. 

*^-1/3  ^ 

Here  E = (1.5)r  , t may  be  any  nonnegative  integer.  It  can  be  shown  that 

[2  sin(X/2)  ]"^F(X)  = (1/12)  (X-rX^)  - (ll/144)rX‘*  - (43/1440)X^  -(l/288)r^X^  +L^(r,X). 
Here  |Lj^(r,X)|  < (0.012)X^  + (0.06)rX^  + (0.012)r^x’  + (0.0012)r^X®  . But 

-ax  n ...n+l  . , ..._,n+l,  ,, , 2 .^2,n+l 

) e X cos(bx)dx  = ri  [ (a-ib)  + (a+ib)  )/2(a  +b  ) 

0 

Hence , 

(4.44)  r e"’^^cos(sX)  (1/12)  (X-rX^)  = - (1/12)  (r^+s^) + (2/3)  r^s^  (r^+s^) = (l/24)r'^, 
0 

when  r = s . 

We  also  have 

(4.45)  (1/72)/"  lllrX^  + (43/10)X^  + (l/2)r^X^  + 144 1 (r , X)  | le‘’^^dX  < 5r‘^  + 308r"^. 
Since  | sinhp | ^ < X ^ for  X < 1 , we  have 

I ;^dX|  < + e-’^‘^)dX. 

e c 

But  e~^  < , 0 < X < e.  Hence, 


(4.46) 


I f"  .1  d.'  1 < 2 (f 


-1.  Sr 


2/1 


+ P 


-15r 


2/3 


-1.35r 


!/3, 


(2/3)  lo<4  r [e  + e 

By  combining  (4.40),  (4.41),  (4 . 44)  - (4 . 46)  and  using  tile  technique  in  the  proof  of  theorem  4.1  to 
estimate  the  remaining  remainder  terms,  we  see  that  Theorem  4.6  holds. 

Theorem  4.6  also  provides  a means  of  generating  G(rh,sh)  by  means  of  a fast  Poisson 
solver  on  a rectangle  using  the  first  three  terms  on  the  right  hand  side  of  Equation  (4.42)  as 
approximations  for  distant  values  of  the  Dirichlet  data. 


We  new  describe  an  efficient  method  developed  in  [24]  of  computing  Gv  for  any  vector 

V defined  on  a square  mesh  S o fi  , with  boundary  mesh  Let  U and  u,  denote  the 

n n h S 3 S 

extension  operators  from  and  3 respectively  to  all  mesh  points  that  are  defined  the 

T 

same  way  as  U . We  are  actually  computing  Ug  G v.  We  first  solve  the  system  of  equations 

T 

U^BU  (J)  = V on  S 

s s 

4»  * 0 on  S 

for  the  potential  • We  then  extend  by  zero  to  all  mesh  points.  We  represent  (p  as 
(4.47)  (J)  = GUgV  + GUggP 

where  p is  an  unknown  vector  defined  on  the  mesh  points  on  9S  to  be  determined.  It  is  easy 
to  see  that 

p = U^gBO.  . 

T T 

The‘  vector  0^  G UggP  can  easily  be  computed  by  one  fast  Poisson  solver  on  S with 

as  the  Dirichlet  data  on  ^S.  Because  of  the  sparsity  of  the  vector  , the  Dirichlet  data 

2 T 

can  be  computed  at  a cost  of  constant  N . ^g''  then  computed  from  (4.47). 

T 

We  now  describe  a method  of  computing  all  three  of  the  vectors  Gv,  Gvjj  and  W Gv  using 

T 

only  two  calls  of  fast  Poisson  solvers.  This  may  appear  to  be  impossible  since  W Gv,  the 
right  hand  side  of  the  capacitance  matrix  equation,  must  be  determined  first  and  the  computation 
of  Gv  alone  requires  two  calls  of  fast  Poisson  solvers.  We  can,  however,  first  compute  and 
store  the  vectors  and  p in  equation  (4.47).  Clearly, 
w'^Gv  - - w'^GUggp. 

T 

But  we  need  only  to  compute  those  mesh  points  that  W Gv  is  defined.  Hence,  the 

T 2 

computation  of  W GUggO  requires  only  constant  n operations.  On  the  other  hand,  ^'^ggP 

T 

GUp  can  be  computed  simultaneously  with  one  call  of  fast  Poisson  solver  with  Ug,G(0ggP+  Up)  as 
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A 


the  Dirichlet  data  on  S . Hence  our  algorithm  of  computing  all  three  of  the  vectors  w'^Gv, 

2 

GUu  and  Gv  requires  only  calls  of  fast  Poisson  solvers  plus  constant  times  n operations. 

In  the  methods  described  so  far,  a fast  solver  on  2N  >■  2N  mesh  points  is  needed  to 

generate  the  discrete  Green's  function  or  its  undivided  differences  on  a N x N mesh  S,  . 

h 

An  alternative  method  is  to  first  generate  the  undivided  differences  of  G on  a N/2  x N/2 

mesh  using  a fast  solver  on  N x n mesh  points.  The  values  on  the  remaining  mesh  points  of 

are  computed  by  using  (4.11)  and  (4.25).  An  accuracy  of  eight  decimal  digits  is  guaranteed 

by  Theorems  4.1  eind  4.2  if  N ^ 60.  A somewhat  less  accurate  but  easier  to  program  method  is 

to  generate  the  values  of  G on  a N/2  x n/2  mesh  and  compute  the  values  on  the  rest  of  S. 

h 

by  using  (4.42)  and  (4.43).  An  accuracy  of  five  decimal  digits  is  guaranteed  by  Theorem  4.6 
if  N ^ 60. 

We  s)iall  assume  in  the  next  two  sections  that  the  G used  in  the  capacitance  matrix 
equation  is  the  discrete  Green's  function  on  the  entire  plane.  The  main  results  in  sections 
5 and  6,  however,  will  also  hold  if  G = is  used  in  equations  (3.2)  and  (3.5).  See 

section  5 of  [31]  for  a discussion  in  this  respect. 
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5 . Spectral  bounds  of  . 

We  shall  show  in  this  section  that  B.  is  uniformly  well  conditioned  in  the  spectral 

h 

norm  as  h ^ 0.  The  following  well  known  lemma  is  crucial  to  the  proof  of  our  main  result. 
Lemma  5,1.  Let  the  symmetric  part  of  a matrix  A satisfy 

(A+a'^)/2  > 6i,  6 > 0 . 


Then 


a'^a  > 6^1. 


Theorem  5.1.  (0.25)  1 < B,  B,  < (7.29)1 

^ h n — 

(0.25)  I < B,  < (13.7)1 

— n n — 

(0.04)  I < B,  < (5.57)1 

— h n — 


for  scheme  I. a 
for  scheme  I.b 
for  scheme  II 


for  all  sufficiently  small  h > 0. 

Proof.  We  shall  first  prove  that  the  following  holds  for  scheme  Ib. 


(5.1) 


^ 1 ®h  «h  • 


(5.2) 


Let  B H B + B,  . We  shall  show  that 
s h h 

nin  {b  (P,P)  - I 1b  (P,Q)  |}  > 1 


Pe3n, 


Qe3fij^,Q;<P 


so  that  (5.1)  holds  because  of  a well  known  Gerschgorin  theorem.  The  inequality 

T 

' i ®h  "h 

will  then  follow  from  Lemma  5.1. 

Let  P e 3n.  . Assume  that  the  local  orientation  of  the  boundary  near  P is  such  that 

h 

for  any  point  P'  £ in  that  neighbourhood,  either  W'  and  N' , the  western  and  northern 

h 

neighbours  of  P',  are  both  in  or  W'  alone  is  in  (Cf))^^.  Let  ttp,  be  the  angle 

£ V4  that  the  normal  through  P'  makes  with  the  x^  axis  in  the  east-west  direction.  By  (3.5) 
and  (3.9),  we  have,  for  P ^ Q , 

(5.3)  B^(P,0)  = 2IG(Wp:W^)  - G(Wp;Q)l  + 2 tanttp [G (Wp; NW^)  - G(Wp;W^)) 

if  P has  only  one  neighbour  in  (Cll)j^.  Here  Yp  denotes  the  immediate  neighbour  on  the  mesh 
for  any  point  P in  the  Y direction.  Similar  expressions  to  (5.3)  are  easily  obtained  when 
Np  is  also  in  when  P = Q . If  P 5 (jh,kh),  Q B (mh,nh)  , then  because  of  transla- 

tional invariance,  G(P;Q)  7 G( | j-m(h, |k-n|h)  . 


Assume  that  0 < a < ti/4  and  that  (912,  ),  , which  denotes  a subset  of  30.  that 

P h loc  h 

contains  a /h  neighbourhood  of  P , can  be  partitioned  into  blocks  as  follows.  Let 
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Then, 


Ig  = {(0,h),...,(0,M^h)} 

= {(kh,M^h-'h),...,(kh,M^^^h)},  k = 1 , 

(-kh,-M_j^h+h)},  k = 1 K^. 


on. ), 

h loc 


K=-K, 


Pel.. 


Note  that  ^ ^ 1 gives  the  number  of  points  in  k = 1, . . . ,K^  while  M - M 


gives  the  number  of  points  in  I , k = 1 , . . . : Mg  = 0. 

Let  denote  the  point  with  x^-coordinate  jh.  From  (5.3)  and  Theorem  4.3,  it  is 

easily  verified  that  \ |b  (P,Q)|  will  remain  essentially  unchanged  for  sufficiently 

smooth  3n  if  tan  is  replaced  throughout  by  tanUp.  Let  a 5 tana^.  Let  P = P^  and 

G(i,j)  = G(ih,jh).  We  shall  assume  that  P ^ P unless  otherwise  stated.  We  easily  verify 

”l 

that  for  P.  e I„, 

3 0 


(5.4) 

(5.5) 

(5.6) 


B (P,P)  = 3+a, 
s 

B (P,P,)  = -4(l+a)  G (0,(i-j|), 
s 3 ^ 

B (P,P„  ) = -2(l+a)  G (0,M  -i)  + (1-a)  G (0,M  -i). 
s X 1 yx  1 


j ^ i,  j / M^  , 


By  (4.39),  we  see  that  for  P^  € Ig  , 

(5.7)  y B (P,P,)  - -2(l+a)[2  G (0,0)  - G (0,M-i)  - g (0,i-l)) 

s j y y 1 y 


+ 2(l+a)  G (0,M  -i)  + (1-a)  G (0,M  -i). 
X 1 yx  1 


For  P.  € I , k * , we  have 

J K 1 

(5.8) 

(5.9) 

(5.10) 


B^(P,P^)  » 2(l+a)  G^(k,j-i-l),  j 

B (P,P.)  = (1+a)  G , (k,j-i-l)  + (1-a)  G^  (k,j-i),  j = M,  . 
s J yy  yx  k+i 

B^(P,Pj)  = -2(l+a)CGy(k^M^-i)  - Gy(k,M^^^-i)I  - (l+a)Gyy {k,M^^^-i-l) 


3 k 


t (1-a)  G^^(k,M^^^-i)  . 


By  Theorem  4.4,  te  see  tliat  each  B (P,P,),  j ^ i,  j < 1,  is  negative.  Hence 

Ki  S3 

(5.11)  I I B_(P,PJ 

k=l  P.  I, 


j k 


j)  ” 2(lta)(Gy(l,M^-i)  t J 

^1  K 


< 2(l+a)  G (1,M  -i) 

y 1 


Similarly, 


(5.12)  I I |b  (P,P.)(  < 2(l+a)  G (l,i-l)  + (1-a)  G {0,i-l). 

k-1  P.el  . ® 3 y y 

j -k 
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By  combining  (5.7),  (5.11)  and  (5.12),  we  eee  tliat 


(5.13)  y Ib  (P,P  )|  < 4(l+a)G  (0,0)  + (1-a)  G (0,i-l)  < (5+3a  G (0,0)  = 5/4  + 3/4a  . 

s I y y y 

Here  we  have  used  (4.5),  (4.6)  and  Theorem  4.4.  Hence,  by  (5.4)  and  (5.13),  we  see  that 

(5.2)  holds  for  P e P P„  . 

U 

The  proof  for  P = P is  quite  similar  and  is  sketched  as  follows. 

M, 


(5.14) 

B (P,P)  = 2+2(l-a) 

G (0,0)  . 

s 

xy 

(5.15) 

y Ib  (p,p.) 

1 = (1+a)  fG  (0,0)  - 

G (0,M  -1) ] + (1-a) [G  (0, 

1)  - G (0,M. )] ; 

y 

y 1 X 

X 1 

(5.16) 

I I 

k=l  P.el,  J 

) 1 = (1-a) [G^(0,0) 

^ J,  ^''xx<’'-^'Vr”i>  - 

k=l 

+ (l+a)[G  (1,0)  + y {G  ) + G 0c,M,  ^ -M  -1) } ] ; 

y yx  k+1  1 yy  Tc+1  1 

(5.17)  I I [b  (P,P  )|  = (lta)[G„(l,M  ) + I {g  ()c,M  +M.  ) + G ^ ()c,M, +M  . -1)  }) 

k=l  P €l  ° y y^  i K yy  i -K 

j -k 

t (l-a)[G^(l,M^)  t - Gy^()t-l,M^tM_(^_^,-l))]. 

By  (4.38)  and  Theorem  4.4,  we  see  that 

(5.18)  y y |b  (P,P.)1  < (1-a)  G (0,0)  + (1+a)  G (1,0). 

k=l  P.el,  ® 5 * y 

3 I' 

Similarly,  by  using  the  identity  (5.20),  we  have 

(5.19)  y y |b  (P,P.)|  < (l+a)G  (1,M,)  + (1-a) (G  (1,M,)  - G (0,M  -1)1  . 

Ic=l  P.el  / ^ ^ ^ ^ ^ ^ 

3 -k 

(5.20)  G ()t,M,+M,-l)  5 G ()c+l  ,M, +M  , -1)  +G  ()c,M,+M.-l)  -G  ()t,M,+M,)  . 

yy  1 -k  yy  1 -k  yx  1 -k  yx  1 -k 

Hence, 

(5.21)  y |b  (P,P.)1  < 2G  (0,0)  + 2G  (0,1)  + aG  (1,M,)  . 

s ] ' — X X y 1 

By  combining  (5.14)  eind  (5.21),  we  have 

(5.22)  B (P,P)  - y |b  (P,P^)|  > 2 + 2(l-a)G  (0,0)  - 2G  (0,0)  - 2G  (0,1)  - aG  (1,M,). 

s s j ' - xy  X X y 1 

Clearly,  the  right  hand  side  of  (5.22)  attains  its  maximum  at  a = 0.  Hence,  (5.2)  holds  for 

P 5 P . The  proof  for  other  choices  of  P is  similar  and  will  not  be  repeated.  We  note 
"l 

that  we  have  assumed  that  B^(P,Q)  ft  0 for  any  Q ^^^^h^loc"  assumption  will  not 

affect  our  estimate  (5.2)  because  each  B (P,Q)  is  either  zero  or  negative  for  P Q. 

s 

Finally  we  remark  that  the  schemes  la  and  Ib  descriljed  in  this  work  are  essentially 

dual  to  the  schemes  I.N.a  and  I.N.b  described  in  (31)  in  the  following  sense.  If  we 

maintain  that  a = tano  does  not  change  its  value  for  the  entire  row  or  column  of  B. 

^ h 
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correspondinq  to  P,  then  B^(P,Q)  is  the  same  for  both  schemes  la  and  I.N.a.  or  for  both 
schemes  Ib  and  I.N.b.  Therefore  we  refer  the  reader  to  | 31)  for  the  proof  of  Theorem  5.1 
for  scheme  I. a.  We  now  proceed  to  prove  the  following  inequality  for  scheme  II. 


(5.23) 


min  {b  (P,P)  - I |b  (P,Q)|)  > 0.41  . 


Ql'P 

Let  P £ . We  assume  the  same  local  configuration  of  irregular  mesh  points  near  P 

h 

as  before.  By  (3.5)  and  (3.7),  we  have  for  ^2  i ^ ^ ^ Q> 


(5.24) 
where 

(5.25) 


Bj^(P,Q)  = -G^(Wp;W^)  - ej  G^(Ep;Wg)  t alG^IWp.-W^)  t e,  G^(Ep;W^)] 


= (1-d^  )/(l+d^ ) , a 5 tan 


Here  G (*;•)  and  G (•;•)  denote  the  forward  undivided  differences  in  the  x,  and  direc- 

X y 12 

tion  with  respect  to  the  second  variable  of  any  mesh  function  G respectively.  A similar  expres- 


wion  involving  e^  which  is  similarly  defined  or  a constant  one  should  be  added  respectively  to 
the  right  hand  side  if  d^  < 1 or  if  P = Q.  Let  and  e^j  denote  the  corresoondino  e^ 

and  e^  respectively. 

Let  P = P..  We  first  estimate  T {|b,  (P,Pj|  +|  B^'cP.P . ) I } . 

" li-j|>3  3 ^ 

* 

Let  B,  (P,Q)  be  defined  by  Equation  (5.24)  with  G replaced  by  its  continuous  analog  which  we 
h 

* •p* 

shall  denote  by  G . B^^  (PtQ)  is  similarly  defined.  By  (5.24)  and  (5.25),  we  have  for  d^  > 1 


and  P ^ 0 

(5.26) 

where 

(5.27) 


(5.28) 


Bj^(P,Q)  = (2/(l+d^)  1 :-G  (W  ;Q)  + (l-a)G  (W  ;W^)  + a G (W  .-NW^)]  + 

= [2/1+dj^))  [G*(W*;W')  - G*(W*;Q)]  + Rq  > 

n-1  k k 

13  = [a/(l+d,)]-[  I (hVDKd.-l)  + (d,+l)'') 

WE  1 k=2  ^ ^ 

•[0/3x^)’‘{-G*(*!Q)  + (l-a)G*(":W^)  + aG*  (•  ;NW^) }]  (W*)  ]} 

+ nth  order  remainder  term. 

Rg  5 (h^a(l-s)/(l+d^)){  atO/SXj)^  G*  (W* ; • ) ) (NW^)  + (1-a)  [ O/ax^)  V (W*;  • ) ) (NW^  ) } . 


Here  W*  is  the  point  were  the  Dirichlet  data  0^  is  given.  W'  is  a point  on  the  mesh  line 
connecting  W^  and  NW^  and  at  a distance  ah  from  W^.  NW^  and  are  respectively  points 

that  can  be  anywhere  on  the  mesh  lines  between  W'  and  W^  and  between  W'  and  We 

shall  assume  that  ^2Q’  analogue  of  d^  for  the  point  0 , is  also  greater  than  1 . In 
that  case,  we  have  a similar  expression  for  bF* (P,Q)  as  that  in  (5.26). 
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We  now  proceed  to  estimate  (PfQ)  ■ It  is  easily  verified  that  if 

P,Q  are  two  points  on  3S2  with  d(P,Q)  < h^,  0 < Y ^ 1 and  ta  is  the  tangent  at  Q to  9f?, 
then 

d(P,ta)  < (K  +o(l))h^^. 

— max 

Here  K is  the  maximum  cibsolute  value  of  the  curvatures  of  9Q.  Hence,  without  loss  of 
max 

* * 

generality,  we  may  assume  that  W lies  on  the  tangent  to  3f?  through  Q , the  point  where  the 

* 

normal  through  Q intersects  with  9Q  and  vice  versa.  Let  r denote  d(W  ,Q)  and  r'  denote 

* *2 
d (W  ,W*).  It  is  easily  seen  that  r'  ^r  if  d(Q  ,Q)  <_h(l+a  )/2.  We  separate  our  discussion 

* 2 

into  four  cases.  The  first  case  is  when  r'  > r and  d(P,P  ) < h(l+a  )/2.  The  maximum  of 

* 2 2 2 2 

log  (rVr)  then  occurs  when  Q coincides  with  Q . In  this  case,  r'  - r « (1+a  )h  . Hence, 
if  P=P^,  Q = Pj,is^j»  then 

0<G*(W*;W')  - G*(W*;Q)  £ (1/4TI)  1 j-i  ( . 

By  (5.28), 

-(1/811)  li-jl'^  1 Rq  1 0 ! if  j > i 

-(1/811)  li-j-l]”^  - *^0  - Ij-i-ll"^  : if  i > j. 

Similarly,  by  (5.27), 

'v'  - '2/'”  li-il"’  • 

Hence,  for  lj-i|  ^ 3, 

(5.29)  (1+d^)  |b*(P,Q)  1 < (1/211)  I j-il"^  + (1/511)  I j-i-ll"^  + (2/11)  I j-ir"*,  if  i>j; 

< (1/211)  Ij-il"^  + (2/11)  I j-il"’,  if  j>i  . 

*2  T* 

Since  d(P,P  ) < h(l’  )/2,  the  estimate  for  (P»Q)  is  the  same  as  that  given  in 

(5.29)  . Hence, 

(5.30)  |B|^*(P,Q)|  + lB*(P,e)l  < (l/H)  1 j-il"^  + (4/11)  1 j-il"'^  + (1/511)  I j-i-ll”^  . 

* 2 * 

The  second  case  is  when  r £ r'  and  d(P,P  ) < h(l+a  )/2.  Let  d^^^  denote  d(Q  ,Q)/h. 
Then  dj^^  > 1/2.  Hence  by  (5.29), 

|b^*(P,Q)|  £ (1/311)  I j-il"^  + (2/1511)1  j-i-ll'^  + (4/311)  I j-i|“‘’  , if  j > i ; 

£ (1/311)  I j-il*^  + (4/311)  Ij-il"*  , if  i > j . 

* 

On  the  other  hand,  the  maximum  of  log  r/r*  occurs  when  Q coincides  with  W^.  In  this 
2 2 2 2 

case,  r -r  * (1-a  )h  . Hence, 

^ ^ ^ ^ 0 0 o 

G (W  ;W’)  - G (W  ;Q)  < (1/4-n)  1 (1-a  )/(l+a  ) 1 I j-i|"  , if  j > i ; 

£ (l/4ii)((l-a^)/(l+a^)]|  j-i-a|^,  if  i>  j ; 
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k 


cuid 


(5.31)  |b*(P,C))|  < [(l/27r)  + a/87T)]  |;j-i|"2  + (2/Tr)  | j-i  | if  j > i ; 

< (1/2  + (l/STT)  I j-i-l|*^  + (211)  I j-il"'’,  if  i > j . 

Tlierefore , 

(5.32)  |b^  (P,Q)  I + |B^(P,e)|  < [(1/3TT)  + (1/211)  + (l/BH)  ] | j-i  + (4/311)  ] | j-i  T'* 

+ (2/1511)  I j-i-lT^,  if  j > i : 

£ [(1/211)  + (1/311))  I j-il"^  + (1/811)  Ij-i-lT^ 

+ 1(4/311)  + (2/11)  ] I j-il"*,  if  i > j . 

The  third  case  is  when  r < r'  and  d(P,P*)  > h(l+a^)/2.  The  estimate  for  |b*(P,Q)|  + 

I T*  1 

+|B^  (P,Q) I is  the  same  as  that  for  the  second  case. 

The  fourth  case  is  when  r > r'  and  d(P,p  ) > h(l+a  )/2.  Both  d^  and  d^^  are  not 

less  than  1/2.  Therefore  by  (5.31)  and  the  above  observation, 

(5.33)  |bJ^*(P,Q)|  + |b*(P,Q)|  < (2/311)  I j-il"^  + (1/1211)  ] j-i  1"^ 

+ (1/1211)  I j-il"^  + (8/311)  I j-il"'*  . 

By  comparing  (5.30),  (5.32)  and  (5.33),  we  see  that 

max  |b^*(p,Q)|  + |b*(P,q|  < (l/ll)|j-ir^  + (4/11)  I j-il"*  + (1/511)  | j-i-1  if  j > i 

< (l/ll)|j-ir^  + (1/811)  I j-il"^  + (1/511)  I j-i-ll"^  + (4/11)  I j-if'' 

if  i > j . 


Hence,  if  d_,  d__  > 1,  then 
2 2Q 

I |bJ[*(P,P.|  + |b*(P,P.)1  < (l/H)  I {2A^  + 8A‘*  + 2/Sk^  + l/8)c^} 
|j-i|>3  " ^ ^ 

< 0.328  . 


Similarly,  it  can  be  shown  that  if  both  d^  and  d^^  are  not  greater  than  1,  then 
(5.34)  y |bJ*(P,P.)|  + |b*(P,P.)|  < 0.677. 

By  Theorems  4. 1-4. 3 and  Table  I on  p.  41  of  [30], 

h *•  3 


I {K(P,P  ) - B*(p,P  )|  + |bJ^(P,P  )|}  < 0.04  . 
i-j|<3  ” J n 3 n 3 


Hence 


(5.35)  I {|b  (P,P  )|  + |b3(P,P.)|)  < 0.717  . 

|j-i|l3  '’3  n 3 

It  remains  to  estimate  B^(P,P)  - y|B^(P,P^)|,  |i-3|  £ 2,  i / 3 . Without  loss  of  generality, 
we  may  assume  that  both  d and  d are  less  than  1.  We  shall  assume  that  PSP.  5 P„  with 

12  i Mj 

^1  > 3 and  i case  when  i / can  be  treated  in  a similar  manner.  We  have 
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B (P,P)  = 2 + 2(l-a)  G (0,0)  + (a)  , 

s xy  0 


g„(a)  H 2 e,  tG  (1,0)  + G (2,0)]  + 2 e,,[G  (0,1)  + a G (1,1)]. 
0 lx  y 2 X y 


For  any  P^  e 1^,  j ^ we  have 

BslP/Pj)  = gj^(a,j)  + 02(3/3)  / 


g (a,j)  = G (0,i-j)  + a G (0,i-j-l)  - a G (l,i-j)  - G (0,i-j)  , 

1 xy  yy  y x 

02(3, j)  = e^^  [G^(l,i-j)  - a G^(2,i-j-l)]  + e^^  (G^  (1 , i- j ) + a Gy(2,i-j)J 


+ e (G  (0,i-j-l)  - a G (l,i-j-2)]  . 
2 X y 


For  any  P.  e I, , we  have 
3 1 

B^(P,Pj)  = g2(a/j)  + g^(a,j)  / 


g,(a,j)  = -G  (0,j-i-l)  + a G (l,j-i-l)  + a G (0,j-i-l) 

3 ■'  X ■'  yy  y 

g^(a,j)  E e^[G^(0,j-i)  + a Gy(l,j-i)]  + e^ l-Gx (0 , j-i+1)  + a Gy(0,j-i+l)] 

+ e.  .(G  (2,j-i)  - a G (3,j~i-l)]. 

By  Theorem  4.4,  gj^(a,j)  is  negative.  It  is  easily  verified  that  02  (a, j)  is  nonnega- 
tive for  0 _<  i-j  £ 2;  gj(a,j)  is  negative  for  a ^ 1/2;  and  g^(a,j)  is  positive.  Moreover, 
for  a £ 1/2, 

2 2 

(5.37)  I Ig  (a,j)  +g2(a,j)l  + I |g2{a,j)  +g^(a,j)| 

i-j=l  j-i=l 

2 2 

' %<^)  I lgj(a/j)l  + I 103(3, j)|. 

i-j=l  j-i=l 

Hence,  for  a < 1/3,  ^3,  M2  ^ have  from  (5.37)  that  the  following  holds. 

(5.38)  B^(P,P)  - I{|b^(P,PJ|  , i / j,  |i-jl  £ 2} 

2 2 

> 2+2(l-a)G  (0,0)  + I g, (a,j)  + I g (a,j) 

’‘y  i-j.l  j-i=l 

> 2+2(l-a)G  (0,0)  - G (0,0)  - 3G  (0,1)  - G (0,2)  + G (0,3) 

- xy  X X X X 


By  considering  all  possible  configurations  of  , i-j  £ 2 , a £ 1/3,  it  can  be  shown 

that  the  constant  1.17  is  always  majorized  by  the  left  hand  side  of  (5.38).  it  is  easy  to  see 

that (5.38)  also  holds  when  P.  / P„  . Hence  by  (5.35),  we  see  that  (5.23)  holds. 

1 M^ 

This  estciblished  the  lower  spectral  bounds  of  for  all  schemes.  To  complete  the 
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proof  of  the  theorem,  we  note  that  the  spectral  norm  of  majorizeri  by  (1/2)  + 

2 

II  Bj^ll^)  . It  therefore  suffices  to  prove  t'.  j following  tvra  inequalities 

(5.39)  II  Bj^ll^  + II  Bj^ll„  ^ if  scheme  I.b  is  used 

(5.40)  II  Ballot,  + II  ®j,lloo  £ 4.72  if  scheme  II  is  used. 

We  first  prove  that  (5.39)  holds.  Without  loss  of  generality  we  may  assume  that 
P e has  only  one  western  neiglibour  W in  (Ci2)j^.  Let  P = (0,0)  and  P^  = (x,y)  , be  in 

3Qj^.  Then  x = ay  + b,  |b|  (l+o(l))h  if  d(P,Pj^)  < /h.  And 

B*(P,Pj)  = -(l/2Ti)log(  (x+h)^+y^]  + (a/2'ti)  log (x^+ (y+h)  + ( (l-a)/2ii)  log (x^+y^) 

= (1/271)  (2h(x-ay)  + (a-l)h^)/r^  + R, 

where  r = d(P,P^)  and  |r|  < (1/277)  [ (2x+h)  ^ + (2y+h)  ^lh^/2r^.  It  is  easily  verified  that 

(5.41)  I |b*(P,P.)|  < 1.25  if  |y|  > 2h 

lyp.  ^ 

3 

By  Theorems  4.1-4. 3 and  the  Table  I on  p.  41  of  (30), 

(5.42)  B^{p,P)  = 2 - (1/2) (1-a)  ; 

h 

(5.43)  y |b^(P,P.)|  < 0.28  + 0.55a/2  ; 

|y|<2h  ^ - 

(5.44)  y |b,  (P,P.)  - B*(P,P,)|  <0.06  . 

|y|i2h  ” ^ h 3 ' 

By  (5.41)-(5.44)  , 

I |Bj^(P,P.)|  < 3.7. 

T 

It  is  easily  seen  that  the  above  inequality  also  holds  when  B^^  is  replaced  by  Bj^.  We  have 
therefore  completed  the  proof  of  (5.39)  . 

Let  P = P^.  By  (5.36) , 

(5.45)  |f.^(P,P)|  < 2+2(l-a)G^y(0,0)  +9^(3)  . 

By  (5.37)  and  (5.45), 

(5.46)  I |b  (P,P  ) |+1b^(P,P.) I < 2+2g  (a)  + 2G  (0,0)  + 6G  (0,1) 

|i-j|<2  1 3 h ] 0 X 

1 4. 

By  (5.35)  and  (5.46),  we  see  that  (5.40)  holds. 
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i 


I 


(5.36) 


where 


where 


where 


B (P,P)  = 2 + 2(l-a)  G (0,0)  + (a)  , 

s xy  0 


g„(a)  H 2 e IG  (1,0)  + G (2,0)]  + 2 e.,]G  (0,1)  + a G (1,1)1. 
0 lx  y 2 X y 


For  any  P^  e 1^,  j ^ M^,  we  have 

B^(P,PJ  = gj^(a,j)  + , 


gj^(a,j)  = G^(0,i-j)  + a Gy^(0,i-j-l)  - a Gy(l,i-j)  - G^(0,i-j)  , 
g2(a,j)  = ej^[G^(l,i-j)  - a Gy(2,i-j-l)]  + e^^  [G^  (1 , i- j ) + a 0^(2, i-j)) 


+ e [G  (0,i-j-l)  - a G (l,i-j-2)l  . 
2 X y 


For  any  P^  ^ ^1'  have 

Bg(P,Pj)  = gj(a,j)  + g^(a,j)  , 


gj(a,j)  5 -G^(0,j-i-l)  + a G^^(l,j-i-l)  + a G^(0,j-i-l) 

g^(a,j)  E (G^(0,  j-i)  + a Gy(l,j-i)]  + e2l-Gx(0,  j-i+1)  + a G^(0,j-i+l)] 

+ e.  .[G  (2,j-i)  - a G (3,j-i-l)]. 

I]  X y 

By  Theorem  4.4,  gj^(a,j)  is  negative.  It  is  easily  verified  that  g2(a,j)  is  nonnega- 
tive for  0 ^ i-j  £ 2;  is  negative  for  a ^ 1/2;  and  g^(a,j)  is  positive.  Moreover, 

for  a £ 1/2, 

(5.37) 


2 2 

I |g^(a,j)  +g2(a,j)|  + I |g3(a,j)  +g^(a,j)| 
i-j=l  j-i=l 


2 2 

' gQ(a)  + I Ig^Ca.j)!  + I |g3(a,j)|. 
i-j=l  j-i=l 

Hence,  for  a < 1/3,  ^3,  M^  ^ 5,  we  have  from  (5.37)  that  the  following  holds. 


(5.38) 


B^(P,P)  - If|B^(P,P.)|  , i fl  j,  |i-jl  < 2} 

2 2 

> 2+2(l-a)G  (0,0)  + I g (a,j)  + I g (a,j) 


i-j.i 


j-i-i 


> 2+2(l-a)G  (0,0)  - G (0,0)  - 3G  (0,1)  - G (0,2)  + G (0,3) 

— xy  X X X X 

> 1.17. 

By  considering  all  possible  configurations  of  P^ , i-j  £ 2 , a ^ 1/3,  it  ceui  be  shown 

that  the  constant  1.17  is  always  majorized  by  the  left  hand  side  of  (5.38).  it  is  easy  to  see 

that (5. 38)  also  holds  when  P.  / . Hence  by  (5.35),  we  see  that  (5.23)  holds. 

1 «1 

This  established  the  lower  spectral  bounds  of  B.  for  all  schemes.  To  complete  the 

h 
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6.  Singular  values  of  and  C . 

We  shall  show  tliat  all  except  a few  singular  values  of  C lie  in  the  interval 

[d,  - e,  d,  + E] , E ^ 0,  where  d,  and  d_  are  the  spectral  bounds  of  B,  . This  is  accomplished 
12  12  h 

by  first  proving  that  the  singular  values  of  cluster  around  that  of  a compact  operator  K . 

Our  main  result  then  follows  as  an  immediate  consequence  of  a well  known  result  in  matrix  theory 
which  will  be  stated  below  as  Lemma  6.8.  We  first  need  some  definitions  from  modern  analysis. 

Let  X denote  a Banach  space  throughout  this  section. 

Definition.  A subset  S C x is  sequentially  compact  if  any  sequence  in  S contains  a convergent 
subsequence  with  limit  in  X . 

Definition.  A family  of  operators  K on  X is  collectively  compact  if  the  set  {k  f : Ilf  I1  <1, 
n in 

fex,  m=l,2,...}  is  sequentially  compact  in  X . 

We  shall  first  assume  that  either  scheme  I. a or  scheme  I.b  is  used.  We  start  by  con- 
structing a family  of  operators  from  in  the  same  way  that  is  done  in  Section  5 of 

[32].  For  completeness,  we  briefly  sketch  this  construction  in  the  following.  Define 

K : C(0,1]  -»  C(0,1]  by 
m 

m 

(6.1)  [K  1 f(t)  = y k(t,t.)f  (t.)  , t.  £ [0,11;  f e C[0,1], 

m ^^^33  3 

where 


(6.2) 


k(t,t.)  = K.  (P.,P.)  + [(t-t.)/(t.  ,-t.)][K,  (P.  ,,P.)  - K,  (P.,P.)]  , t.  ,<  t < t,  . 
3 h 1 g 1 1-1  1 h 1-1  g nig  i-l  — — i 

C(0,11  is  the  Banach  space  of  continuous  functions  on  [0,1].  The  t^,  i = l,...,ro 


are  defined  as  follows.  Let  (fi.ilJ  be  a smooth  parametrization  of  35).  Then  (ij)  (t^)  ,i(;(t^) ) is  the 

closest  point  on  30  to  P.  e 30  which  is  on  the  normal  through  P..  When  t is  very  close  to 

1 h 1 

0 or  1 , the  k(t,tj  in  (6.2)  should  be  adjusted  slightly.  See  [32]  for  the  details.  We  can 

construct  by  the  same  procedure  a family  of  operators  {K'}  from  {kJ  }.  Let  K ^ k'K 

in  n s m ID 

T T 

Lemma  6.1.  The  nonzero  eigenvalues  of  3nd  coincides  with  that  of  K^, 

K +K'  and  K respectively, 
mm  s 


Proof : See  e.g.  Lemma  5.2  in  [32]  , 


B 


Lenima  6.2.  Let  P and  Q be  two  points  in  30^^  with  d(P,Q)  = h^,  B £ 1/2.  Let  P and  Q 
the  closest  points  on  30  to  P and  Q respectively.  Then 
lC^(P,e)  = 2[3g  /3v^,)(p  ;Q  )hseca^  + 0(h  ‘^)  . 

Proof ; Essentially  the  same  as  that  of  Lemma  5.4  in  (31). 
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Lomma  6.3. 


The  families  of  operators  (k  and  (k  } are  collectively  compact  on  C(0, IJ. 

mm  s 

Proof : Essentially  the  same  as  that  rf  Lemma  5.5  in  [32]. 

T T 

Lemma  6.4.  K^f  Kf,  K^f  K f an-:!  K^f  -►  K Kf  for  each  f € C[0,1]  where  K is  the  compact 
integral  operator  defined  by 


(Kf)(t_)  = 2 j [3gV3v„]  (P;Q)f  ds„ 
where  P : ((S’ (tp)  .ij;  (tp)  . 

Proof : Essentially  the  same  as  that  of  Lemma  5.6  in  (32]  . 

In  order  for  the  above  theorems  to  apply  in  the  case  when  scheme  II  is  used,  we  scale 
the  matrix  in  that  case  as  follows.  The  rows  of  that  correspond  to  irregular  mesh 

points  that  have  one  or  two  neighbours  in  multiplied  with  (l+d^^)  or 

(1+dj^)  (l+d^)  (l+dj^+d^)  ^ respectively.  It  is  easily  verified  that  Lemmas  6.3  and  6.4  hold  for 
scheme  II  if  is  constructed  from  the  scaled  It  will  be  shown  after  theorem  6.1  that 

such  a scaling  is  not  e-sential  and  our  main  results  will  hold  even  without  it. 

Lemma  6.5.  Let  {k  } be  collectively  compact  on  X : K f Kf  for  each  f e X.  Given  e > 0, 

n n 

let  y^,  with  algebraic  multiplicities  m^,  i = 1,...,N  be  the  eigenvalues  of  K with  absolute 

* ★ 

values  greater  them  or  equal  to  c > 0.  Then  there  exist  positive  numbers  N and  e < e 

* * 


such  that  for  all  n ^ N each  e neighbourhood  of  contains  exactly  m^  eigenvalues  of 

K while  all  the  other  eigenvalues  of  K lie  in  an  e -neighbourhood  of  zero, 
n n 

Proof ; This  is  an  immediate  consequence  of  Theorem  4.8  on  p.  65  of  (1].  See  also  Chapter  4 of 
(301.  By  combining  Lemmas  6.1,  6.3,  6.4  and  6.5,  we  easily  have  the  following. 

Theorem  6.1  Given  e > 0,  there  exists  a positive  integer  N such  that  for  all  h > 0,  all 
except  N singular  values  of  Kj^  lie  in  (0,el. 

Lemma  6.6.  Let  C * AB,  where  A,  B and  C are  arbitrary  matrices  with  singular  values 

a,  > o,  > . . . > a > 0,  6,  > 6-  >. . .>  6 >0  and  Y,  > Y-  Y >0  respectively,  then 

1 — 2 — — m — 1 — 2 — — m — 1 — 2 — — m — 

i “i+l® j+l  ' ^’3  positive  integers. 

Proof:  See  e.g.  Exercise  28  on  p.  89  of  (23]  . An  immediate  consequence  of  Lemma  6.6  is  that 

Theorem  6.1  holds  in  the  case  when  scheme  II  is  used  even  if  the  matrices  C or  are  not 

scaled  by  the  scaling  described  just  before  Lemma  6.5. 

Lemma  6.7.  If  D « A+B,  where  A and  B are  as  in  Lemma  6.6,  and  S > 6_  > . . . > 6 >0  are 

— ■ 1 “ 2 ■“  m ” 

the  singular  values  of  D , then 

"^i+j+l  - “l+l  ®j+l  , l,j  positive  integers. 
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I 


Proof : 


See  e.g.  Exercise  30  on  p.  89  of  (231. 


Theorem  6.2  Let  dj  and  d^  be  the  spectral  bounds  of  Then  given  G > 0,  there  exists 

a positive  integer  N independent  of  h such  that  all  except  N singular  values  of  C lie  in 
[d^-e .d^+t] . 

Proof ; An  immediate  consequence  of  Lemma  6.7  and  Theorem  6.1.  See  also  Theorem  5.3  in  (32).  In 
the  following,  ll  ||  shall  denote  either  the  spectral  norm  of  a matrix  or  the  Euclidean  norm  of  a vector . 

•k 

Lemma  6.8  Let  U be  the  extension  operator  from  52^  u ot  all  mesh  points  that  is  defined 

*T 

the  same  way  as  U . Suppose  that  0 G Vy  ^ 0 for  any  nonzero  m-vector  y defined  on 

Then  C is  nonsingular.  Moreover,  if  ||  U G Vy||  ^ C |(  yjl  /|(  A (j  for  any  m-vector  y 

^ 11 

then  II  C ^11  £ K (Aj^j^)/Cj^,  where  x (A^j^)  is  the  spectral  condition  number  of  A^j^  with  respect  to 

the  norm  ||  || 

T T 

Proof : Let  Au  = v = UU  v with  U v in  the  range  of  C be  the  equation  we  are  solving. 

From  Section  3,  we  see  that  u = GVy  is  a solution  of  Au  = v if  y satisfies  Equation  (3.6). 

Suppose  C is  singular  so  that  there  exist  two  distince  solutions  y^  and  y^  of  Equation  (3.6). 

»T 

Let  y^  = iJj^-y2  • Then  ^^Vy^^  = 0.  Because  of  the  reducible  structure  of  A , GVy^  = 0. 

*T 

This  contradicts  the  assumptions  that  Aj^j^  is  nonsingular  and  U GVy^j  0.  Moreover,  if 
II  U*’'GVy||  > C^ll  y||  /ll  A^^ll  , then 

II  Aiill  II  u\|l  > ||u*’'GVy||  > c^ll  y||  /|1  a^^||  . 

The  lemma  easily  follows. 

Definition.  A scheme  of  interpolating  boundary  conditions  is  said  to  be  admissible  if  its  cor- 
responding coefficient  matrix  Aj^^^  of  the  discrete  problem  is  nonsingular  and  k(A^^)  ^ 

-2 

constant  h 

Lemma  6.9  Let  C*  and  A*  denote  respectively  the  capacitance  matrix  and  the  coefficient 
matrix  of  the  discrete  problem  for  a certain  scheme  of  interpolating  boundary  conditions.  Suppose 
that  both  C*  and  A*  are  nonsingular.  Then  C is  nonsingular  for  any  admissible  scheme  of 
interpolating  boundary  conditions.  Moreover,  if  ||  c ^11  1 <=2  \l^^  ' 

then  II  C ^11  1 (Aj^^)  . 


Proof:  '^e  first  claim  that  if  both  C and  A 


are  nonsingular,  then  there  exists  no  u ^ 0 


such  that  U GVy 


Suppose  this  is  not  so,  then  there  exists  a p / 0 such  that 


♦ Y * * 

U GVp  - 0.  Let  V = UU  V where  C p s=  U v . Since  C is  nonsingular,  U v ^ 0.  But 

★T  * *T  * * *'p  * 

U V = U AGVy  = A U GVVI  * 0 . 

This  proves  our  claim.  By  Lemma  6.8,  C is  nonsingular.  Suppose  now  that  ||  C ^ II  £ <=2 


A II  £ c^ll  A^^ll  . Let  Au  = V = UU  V.  Clearly, 

II  Pill  |lc*-ll|  1 l|c‘-^||  II  A* 


u gvpII  . 


Hence, 


II  U*'^GVp||  1 II  U||  /C2C3II  A^J|. 

The  lemma  easily  follows  from  Lemma  6.8. 

Definition  U is  said  to  be  in  3 (U)  if  the  associated  integral  operator  K defined  by 
Equations  (2.2)-(2.3)  is  such  that  K + > -SI. 

Lemma  6.10.  All  ellipses  with  thickness  b/a  > 1/3  are  in  3(1).  Here  a and  b are  re- 
spectively the  major  and  minor  axes  of  the  ellipses. 

Proof : An  immediate  consequence  of  (2.5). 

Theorem  6. 3 Let  G = be  the  discrete  Green’s  function  used  in  equations  (3.2)  and  (3.5). 

Then  the  capacitance  matrix  C is  nonsingular  and  ||  C ^||  1 constant  h”"^  for  some  positive 
integer  q independent  of  h . 

* * 

Proof:  It  suffices  to  find  a pair  <C  ,A  > that  satisfies  the  hypothesis  of  lemma  6.9.  Assume 

that  the  difference  equations  are  already  preordered  in  such  a way  so  that 


where  the  first,  second  and  last  rows  of  B in  block  form  correspond  to  the  coefficient  matrices 
of  the  difference  equations  on  abd  respectively. 

Suppose  that  in  forming  A we  use  a zero  order  interpolation  of  Dirichlet  data  on  30 
at  to  obtain  the  equations  on  3 0j^.  Partition  v,A  and  G in  the  same  way  as  B.  We  obtain 


so  that  the  capacitance  matrix  which  we  now  denote  by  satisfies 
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It  is  easily  seen  that  A,,  is  the  coefficient  matrix  of  a discrete  exterior  Neumann  problem, 

IN 

with  the  normal  derivative  approximated  by  the  first  order  scheme  described  on  p.  203  of  [13). 

Hence  A.,  is  nonsingular  since  both  B,,  and  A,.  are  nonsingular  and  A is  reducible. 

N 11  IN  N 

Clearly, 

(6.5)  A„  = B + Uv"^  - Uu'^B  . 

N 

Let  u , be  the  solution  of 
N 

Vn  = ' 

where  f^  is  any  mesh  function  that  vanishes  outside  Suppose  we  make  the  Ansatz  that 

(6.7)  u,  = B ^ Up 

N 

where  p satisfies 

(6.8)  u'^  A„  b"^  up  = u’’  f . 

N N 

T -1 

By  (6.5),  (6.7)  and  (6.8),  it  is  easily  seen  that  (6.6)  is  satisfied.  Let  C 5 U A B U.  It 

N N 

is  clear  that  C is  nonsingular.  By  (6.4)  and  (6.5) , we  have 
N 


so  that  is  nonsingular.  Moreover,  using  an  argument  similar  to  the  proof  of  lemma  6.8,  we 

have 


for  some  positive  integer  q independent  of  h . The  Theorem  easily  follows. 

Theorem  6.4  Let  SJ  e 2f(l).  Assume  that  the  G in  equation  (3.5)  is  the  discrete  Green's  func- 
tion of  the  entire  plane.  Then  ||  C ^||  ^ constant  h ^ as  h -►  0 for  any  admissible  scheme 
of  interpolating  boundary  conditions.  Moreover,  1|  C ^|1  ^ constant  as  h -►  0 if  either  scheme 
I. a or  scheme  I.b  is  used  or  if  scheme  II  is  used  and  (2  £ 3(0.4). 

Proof;  Let  tJ  c 3(1).  Assume  that  either  scheme  I. a or  scheme  I.b  is  used  for  interpolating 

T 

the  Ixpundary  conditions.  By  (5.1),  Bjj  + i I-  By  assun^Jtion,  there  exists  an  t > 0 such 
T 

that  K + K ^ -I  + E . By  lemmas  5.1  and  6. 3-6. 5,  we  see  that  for  sufficiently  small  h , 


-29- 


^ -I  + e/2.  Hence,  C + c'^  ^ e/2  and  ||  C ^||  £ constant  as  h ->■  0 . Similarly,  it 
can  be  shown  that  ||  c ^||  ^ constant  if  scheme  II  is  used  and  fl  e 3(0,4).  By  lemma  6.9, 

II  c II  ^ constant  • h for  any  admissible  scheme  of  interpolating  boundary  conditions. 


§7 . Convergence  of  conjugate  gradient  iteration. 


T 

Let  b denote  the  right  hand  side  of  the  capacitance  matrix  equation  multiplied  by  C . 

T 

Let  Q denote  C C.  We  are  concerned  with  solving  Qp  = b by  the  conjugate  gradient  method. 
Detailed  exposition  of  the  method  can  be  found  e.g.  in  [11],  [16],  (17),  [18]  and  (26).  A brief 
description  of  the  method  plus  a simple  extension  of  the  known  results  in  the  above  references 
can  be  found  in  Section  6 of  (32).  It  will  be  assumed  that  the  readers  are  familiar  with  the 
results  in  [ 32) . 

Let  denote  the  vectors  approximating  the  solution  P generated  by  the  conjugate 

gradient  process.  Let  R denote  the  set  of  real  numbers  and  L denote  the  set  of  m vectors. 

m 

2 

Let  Z * R -►  R and  E : L^  R be  defined  respectively  by 

Z(a,b)  = {(1  - v^)/(l  + 

E(Pj^)  r (1/2)  (pj^  - P)'^0(P^  - P)  . 

It  is  shown  in  (321  that  the  following  holds. 

T 

Theorem  7.1  Let  K and  K,  be  the  spectral  condition  numbers  of  0 and  B.  B.  respectively. 

— X X\  tx 

* ^ 

Let  d and  d*  denote  the  smallest  and  largest  eigenvalues  of  B,B,  resepctively.  Then  given 

h h 

f.  > 0,  there  exists  a positive  integer  independent  of  k and  h such  that 

E(P|^)/E(Pp)  < min{4Z(k,2k)  , 4Z  (kj^-2e/d ' , 2k-2N) X } - 

Here  “ “a’'  n|l->i/X.l,  i = 1,...,N,  where  X.,  i = X,...,N  are  the  N eigenvalues  of 

d*<X<d’  „ ^ ^ 

Q that  lie  outside  of  (d  -f,d'+e)  . 

Corollary  7.1  Let  o = bo  used  in  equation  (3.5).  The  number  of  iterations 

needed  to  reduce  E(pj^)/E(Pp)  to  a given  accuracy  can  grow  no  faster  than  constant  "log  m as 

h 0 . 

Proof.  By  Theorem  6.3,  |x(^)|  ^ h where  k is  a constant  independent  of  h.  The  corollary 

is  therefore  an  easy  consequence  of  Theorem  7.1. 

* 

Corollary  7.2  Let  r 3 (p) , g = l if  either  scheme  I. a or  I.b  is  used;  g = 0.4  if  sheme 

II  is  used.  Then  the  number  of  iterations  needed  to  reduce  E(p|^)/E(Pg)  to  a given  accuracy 

stays  constant  as  h 0 if  the  G in  (3.5)  is  the  discrete  Green's  function  on  the  entire  plane. 
Proof.  By  Theorem  5.1  and  Theorem  6.4,  C is  uniformly  well  conditioned  in  the 
spectral  norm.  The  corollary  is  therefore  an  immediate  consequence  of  Theorem  7.1. 
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§8  Survey  of  previous  work  on  capacitance  matrix  methods. 

R.  W.  Hockney  in  [201  and  (211  described  a method  of  this  type  which  can  be  used  for  the 
solution  of  the  interior  Dirichlet  problem  for  Laplace's  equation.  His  capacitance  matrices  are 
always  positive  definite  symmetric.  His  method  thus  corresponds  to  a single  layer  Ansatz  for 
the  Dirichlet  problem.  Buzbee,  Dorr,  George  and  Golub  used  a similar  method  in  [8J . They  made 
the  Ansatz 

u=B^v+B^Uwp  , 

when  B is  nonsingular.  Here  W is  a m ^ m nonsingular  metrix.  The  choice  W = I gives  the 
Woodbury  formula. 

Proskurowski  and  Widlund  introduced  the  double  layer  Ansatz  in  [29] . The  algorithm  used 

in  their  work  differs  from  the  one  used  here  only  in  the  discrete  Green's  function  G and  the 
T 

W matrix.  No  theoretical  analysis  was  presented  in  [29].  In  [2]]  the  author  analyzed  the 
method  for  the  Neumann  problem.  The  algorithm  used  in  [32]  is  similar  to  the  one  used  by 
George  in  [15]  which  corresponds  to  solving  the  single  layer  Ansatz  of  the  Dirichlet  problem 
in  an  iterative  imbedding  fashion. 


§9  Numerical  experiments 

The  results  in  this  section  were  obtained  on  the  CDC  7600  at  Lawrence  Berkeley  Laboratory 
The  model  problem  is  the  Laplace  equation  on  ellipses  with  y = (a-b)/(a+b),  where  a,b  are  the 
half  axes.  E(u)  = ||  u-u  where  u is  the  true  solution  of  -Au  = 0 on  ti  , u = 1 on  30  . 

The  mesh  size  h = 1/32.  The  number  of  iterations  of  the  conjugate  gradient  method  is  denoted 
by  n N(r)  denotes  the  normalized  norm  of  the  residuals  which  is  the  L^  norm  of  the  residual 
divided  by  the  square  root  of  points  in  0.  The  numbers  given  for  E(u)  are  acturally  upper 
bounds  that  describe  the  number  of  accurate  digits  only.  The  capacitance  matrix  is  generated 
explicitly  and  the  discrete  Green's  function  on  the  plane  is  used  in  (3.5). 


TABLE  I 


Scheme 

I. a 

Scheme 

I.b 

Scheme 

II 

n 

y 

N(R) 

E(u) 

N(R) 

E(u) 

N(R) 

E(u) 

4 

0.2 

— 

— 

— 

— 

3.9-04 

1.0-03 

5 

0,2 

— 

— 

— 

— 

2.1-04 

1.0-03 

4 

1 

1.5-04 

1.0-03 

8.7-03 

1.0-02 

3.5-04 

1.0-03 

5 

1 

1.0-04 

1.0-03 

4.2-03 

1.0-02 

1.6-04 

1.0-03 

In  Table  I we  see  that  typically  it  takes  four  iterations  to  achieve  three  digits  accuracy. 

2 

The  operation  count  of  the  conjugate  gradient  routine  is  therefore  approximately  64n  . The  total 

operation  count  (not  counting  that  of  setting  up  the  matrix  c)  is  therefore  approximately 

2 2 2 2.. 

5n  logn  + 80n  for  the  Laplace's  equation  and  lOn  log  n + 120n  for  the  Poisson  equation. 
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